Bose glass and Mott insulator phase in the disordered boson Hubbard model 
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We study the Villain representation of the two-dimensional disordered boson Hubbard model via 
Monte Carlo simulations. It is shown that the probability distribution of the local susceptibility has 
a 1/x^-tail in the Bose glass phase. This gives rise to a divergence although particles are completely 
localized here as we prove with the help of the participation ratio. We demonstrate the presence of 
an incompressible Mott lobe within the Bose glass phase and show that a direct Mott-insulator to 
superfluid transition happens at the tip of the lobe. Here we find critical exponents z = 1, ~ 0.7 
and 77 ~ 0.1, which are reminiscent of the pure three-dimensional classical XY model. 



PACS numbers; 67.40-w, 74.70.Mq, 74.20.Mn 

At zero temperature two-dimensional systems of in- 
teracting bosons can show a quantum phase transition 
from an insulating phase to a superconducting phase 
@,|. Such a transition can be observed experimentally in 
granular superconductors |^ and in ''He-films absorbed in 
arogels [Q . By tuning a control parameter like the disor- 
der strength or the chemical potential the bosons become 
localized in a so called Bose glass phase that is insulating 
but gapless and compressible. A huge theoretical effort 
has been undertaken to shed like on the universal prop- 
erties of this superconductor-to-insulator transition. In 
two dimensions the model has eluded successful analyt- 
ical treatment, which necessitates numerical methods as 
quantum Monte Carlo simulations , real space renor- 
malization group calculations pO| or strong coupling ex- 
pansion |Tl| ]. 

Apart from this generic transition the Bose glass phase 
itself has a number of universal features that are rele- 
vant for experiments. Since it is gapless various zero- 
frequency susceptibilities will diverge which is remi- 
niscent of the quantum Griffiths phase occuring in ran- 
dom transverse Ising systems p^-p7|, where a continu- 
ously varying dynamical exponent parametrizes the oc- 
curing singularities. Moreover, for weak disorder a dif- 
ferent transition, directly from a superconducting to a 
Mott-insulating phase might occur This scenario 

emerges also from recent theoretical considerations fl^ 
and would establish a new universality class different 
from the one investigated in 

In this letter we address these two questions in a nu- 
merical approach. We report on results obtained by ex- 
tensive quantum Monte Carlo simulations of the disor- 
dered boson Hubbard model (BH) with short range in- 
teractions in two dimensions, which is defined by 
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where (ij) are nearest neighbor pairs on a square lat- 
tice, af (oi) are boson creation (annihilation) operators. 



Ui — af ai counts the number of bosons at site i, U is 
the strength of an on-site repulsion and /ii is a random 
chemical potential. We are interested in the ground state 
properties (i.e. at temperature T = 0) of (^. By using 
standard manipulations |^] we rewrite the ground state 
energy density of (|l|) as a free energy density of a classical 
model 
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V-Ji,t=0 

In exp(-5{J}), 
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where the integer current variables Jfj, Jf^ and JJ^, live 
on the links of a (2-f l)-dimensional cubic lattice of linear 
size L in the two space directions (with coordinates i = 
{x,y)) and Lt cx 1/T in the (imaginary) time direction 
(with coordinate t). Ultimately one has to perform the 
limit Lr —>■ GO (i.e. T — > 0). The current vector Ji^t = 
iJi,t: Ji,t' '^i,t) has to be divergenceless on each lattice site 
(i,t) as indicated. The classical action S is given by 



(3) 



The coupling constant K acts as a temperature and cor- 
responds to t/U. Note that the mapping from (|l|) to (||) 
involves various approximations and we stress right 
from the beginning that we report exclusively on results 
for the classical model (H). However, as far as univer- 
sal properties are concerned, we expect them to be valid 
also for (|l|). The random part Vi of the local chemical 
potential is distributed uniformly between —A and -l-A. 
All results are disorder averages over at least 500 sam- 
ples, obtained by Monte Carlo simulations of the classical 
model (^) with an appropriate heat bath algorithm at 
classical temperature K. Details of the calculations will 
be pubhshed elsewhere |Q. 

In mean field theory one expects the following K ~ fi 
phase diagram For A < 0.5 there is a superfluid (SF) 
phase at large K, a Bose glass (BG) phase at small K 



1 



and a sequence of Mott-insulator (MI) lobes embedded 
into the BG phase centered around K = 0, ^ integer. 
For A > 0.5 the Mott lobes vanish and only the BG and 
SF phases remain. In this case the SF-BG transition is 
generic everywhere along the phase separation line and 
has been investigated extensively in Q in two dimensions 
at the point K ~ Kc{A = 0.5,^ ~ 0.5). The nature of 
the transition at the tip of the Mott lobes (i.e. at ^ = n 
and K = K'^{A, n) for A < 0.5) is not clear and under 
discussion in the literature |p|,p^,^JTl[| . 

Our first goal is to shed light on the Bose glass phase it- 
self. It has been argued ||] that here the density of states 
at zero energy does not vanish, leading to a divergent 
superfluid susceptibility, although the correlation length 
is finite. On one hand, this behavior is reminiscent of 
the quantum Griffiths phase in random transverse Ising 
systems fl^-IT^. On the other hand we demonstrate in 
this letter that the BG phase is different from the Grif- 
fiths phase in the following respect: whereas in the latter 
strongly coupled clusters lead to a divergence of varying 
strength with varying coupling constant, essentially fully 
localized excitations give rise to a uniform, logarithmic 
divergence in the former. 

We study the local superfluid susceptibility, which is 
defined by Xi — Ylit=i'~'t i^) with the imaginary time 



■ •) mean a ther- 



autocorrelation function Cf{t) = (exp{- 
Jj^j, —^J'^)}) where the angular brackets (■ 

modynamic average. Note that Cf{T) corresponds to the 
local (imaginary time) Greens function (ai(T)a^(O)) in 
the original BH model (|l|) and the local susceptibility is 
simply its (zero frequency) integral. 
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FIG. 1. The probability distribution P(ln Xioc) of the local 
susceptibility for various values of A' (A = 0.5, ^ — 0.5, the 
BG-SF transition is at Kc = 0.247). The system size is L = 6 
and Lt — 200. For K = 0.19, which is deep in the Bose glass 
phase, also data for L = 4 and L = 10 are shown, which is 
indistinguishable from L — 6. 

The probability distribution P(lnx) is shown in fig. |l| 
for the case A — 0.5, from which wc conclude that 



with z = d = 2 throughout the BG phase. We have 
chosen the notation of ref. ||l^,|l^ in order to demonstrate 
that the dynamical exponent z that is characteristic for 
a Griffiths phase in random transverse Ising models 
can also be defined in the present context and is constant 
here. Note that here z ~ d in the BG phase and at 
the critical point, although the two exponents have their 
origin in different physics . We also looked at weaker 
disorder A = 0.2, where MI lobes are present. As soon 
as one enters the latter, the distribution P{x) is chopped 
off at some characteristic value inversely proportional to 
the non- vanishing gap in the MI phase. This implies 
furthermore that the BG phase is indeed gapless ||^. 
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FIG. 2. The participation ratio pl{K) as a function 
of K for various system sizes (A = 0.5, fj, = 0.5). Within the 
Bose glass phase {K < Kc = 0.247) pl approaches a constant. 
The insert shows a scaling plot with i/ — 0.9 and y = 1.0. 

The relation (^) could be obtained by setting the hop- 
ping matrix element t to zero in (|^) , which yields a com- 
pletely local Hamiltonian. This lets us suspect that the 
fact that z does not vary within the BG phase is due to 
the local nature of the low lying excitations. To further 
clarify this point we try to quantify the degree of local- 
ization of the latter. However, since it is not possible to 
obtain these excitations directly in the representation we 
use, we introduce a participation ratio pit] 
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In P(ln x) = In X 



const. 



(4) 



for the spatial density distribution pi — L.^^ ^^{Jli — 
Jl'i ) of an additional particle, i.e. here we work with two 
replicas a and (3 of the system, one with fixed particle 
number N and the other with N + 1 (iV = L^/2 for the 
transition that corresponds to /i = 1/2). [• • -Jav denotes 
a disorder average. One expects pl — C(l) if the extra 
particle is localized and pl — 0{N) if it is delocalized. 
The result is shown in fig. ^ where we see very clearly. 



2 



that the additional particle becomes completely localized 
within the BG phase, most probably at those sites, which 
allow for an extra particle, i.e. « (since we are at 
fi = 1/2). Moreover the insert shows that pl{K) satisfies 
the following scaling relation for fixed aspect ratio Lt/L^ 
at the generic SF-BG transition 



Pl{K) 
L2 



L-yq{SL 
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with S = {K — Kc)/Kc the distance from the critical 
point {Kc = 0.247), = 0.9 ± 0.1 and z = 2 as in |§ and 

y 



1.0 ±0.1. 
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FIG. 3. The average particle number rn and compress- 
ibility K as a function of the chemical potential fi (A — 0.2, 
K = 0.19). At/i = 0.5 one is deep in the Bose glass phase 
and m = 1/2. The lower insert shows the compressibility: for 
large Lt the drop in k becomes sharper. The upper insert 
shows how the plateau (i.e. the Mott lobe) vanishes in the 
vicinity of K = K'^{A = 0.2, = 1). 

Now we consider the MI lobes within the BG phase. 
We choose A — 0.2 and explore the MI-BG boundary by 
varying the chemical potential between 0.5 and 1.5 with 
fixed K < Kc{A = 0.2, /i = 0.5) = 0.20. The latter point 
corresponds to the generic BG-SF transition studied in 
0, but now with weaker disorder. We checked that one 
indeed gets the same critical exponents as for A = 0.5. In 
fig. ^ one sees that with increasing ^ the average particle 
density per site increases monotonically until it saturates 
in a plateau at = 1. The plateau region indicates the 
boundary of the MI phase centered around /i = 1 with 
exactly one particle per site. There is only a weak size 
dependence, at least as long as one is deep inside the BG 
phase, where the correlation length is very small. The 
insert shows the compressibility 
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where N = t) '-^It ^^e total number of particles. 
Obviously the compressibility vanishes as soon as the 
plateau, i.e. the MI phase, is reached. We note that 
we observe extremely strong sample to sample fluctua- 
tions in the compressibility, which necessitated an ex- 
tensive disorder average (5000 samples). With increas- 
ing K the plateau region shrinks until it vanishes at 
K'^{A = 0.2, /i = 1) « 0.325. This indicates the tip 
of the lobe on which we focus now and for which the uni- 
versality class might be different from the generic case. 
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FIG. 4. The stiffness p(0)L at the tip of the lobe (A = 0.2, 
= 1.0). The aspect ratio is constant for z = 1, the insert 
shows a scaling plot with Kc ~ 0.325 and f = 0.7. 

We fix fi — 1 and vary K. Coming from the SF-phase 
we first analyze the finite size scaling behavior of the 
superfluid stiffness 
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with Tlx — l/L J^^ the winding number. Since we 

do not know the dynamical exponent z we hypothesized 
z = 2 (as for the generic case) and z ~ 1 (as in the pure 
(A — 0) case). For both we performed runs with constant 
aspect ration Lr/L^ , and it turned out that for z = 1 the 
best data collapse could be obtained and that only this 
value is also compatible with the correlation function re- 
sults discussed below. In fig. ^ we show our results for 
z — 1, where one gets a clear intersection point of Lp at 
K'^ ~ 0.325 ± 0.002. This value is very close to the corre- 
sponding value for the pure case K'p^'^° = 0.333 ± 0.003 
JlSf and indicates that the tip of the lobe depends very 
weakly on the disorder strength or is even independent of 
it over some range |l^ . In the latter case the critical ex- 
ponent u might escape the inequality u > 2/d |^ at this 
special multicritical point p^ , p3[ |, since then variation of 
the disorder cannot trigger the transition as required in 
p2| . Indeed, the insert of fig. 4 shows a scaling plot which 
yields — 0.7±0.1 < 2/d — 1, which agrees well with the 
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pure value (see below). These results yield a consistent 
picture, nevertheless we should mention that one cannot 
strictly exclude the possibility that our data are not yet 
in the asymptotic scaling regime and the exponent v we 
estimate is only an effective exponent for small length 
scales. 

From the data for the imaginary time correlation func- 
tion Ci(i) and the spatial correlation function C^^fr) — 
[(exp{--^ ELi(l/2 + Jt+xe^,i)\)\av shown in fig. g we 
get firm support for z = 1: The ratio of the decay ex- 
ponents Dx, and y^- for Cx and C+(t) = [C'i^(''')]av, re- 
spectively should equal z and we find DxIVt ~ 1-0(1), 
roughly independent of how we scale h^- with L. From 
yx— d + z — 2 + riwe get 77 = 0.1 ± 0.1. 
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FIG. 5. The correlation functions Cx{r) and C+(r) 
for system size 10x10x10 at the tip of the lobe (-R'c=0.325, 
A=0.2, /i=1.0). The dotted lines are least square fits to 
Cx{r)^a{ry-+{L-rf'^) and C+(r)=a'(r''"-h(Lr-r)''-), which 
gives yx~~l-ll and j/t=-1.08, i.e. z=yx/yT~l-0. The insert 
shows the compressibility k for A=0.2, n—l.O. 

Finally the insert of fig. || shows the compressibility 
at the above studied transition, and we observe that it 
vanishes here, too. In particular the data scale well ac- 
cording to K = L~^k{L^/'^ 5) for systems with L = L^., i.e. 
constant aspect ratio for z = 1. Hence, for weak disorder 
(A < 0.2) we find at the tip of the lobe a direct SF to 
MI transition possibly within the same universality class 
as the pure model, since our estimates for z, rj and v are 
numerically indistinguishable from the values for the the 
pure XY model in 3d, which are z = 1, 77 = 0.033(4) and 
u = 0.669(2) ||. 

As one can see from fig. ^ and the insert of fig. || there 
is no sign of a first order transition at the tip of the lobe, 
as has been suggested in |l^. Moreover, our conclusion 
disagrees with the mean-field prediction ||] of an inter- 
vening BG phase between MI and SF phase at the tip 
of the lobe for weak disorder. For stronger disorder the 
scenario might change: for instance at A = 0.4 we esti- 
mate z « 0.4, which is possibly only an effective exponent 
and the compressibility does not vanish immediately be- 



low the transition from the SF phase. This indicates the 
existence of a threshold value for the disorder strength: 
only above this threshold the mean-field prediction might 
be correct [|o|jl|. 

To conclude, we have shown in this letter that the Bose 
glass phase in the disordered boson Hubbard model and 
the Griffiths phase in random transverse Ising models 
are closely related and that the gapless low energy ex- 
citations are fully localized in the BG phase. Moreover, 
we presented evidence that the transition for commen- 
surate boson densities is directly from a Mott insulating 
phase to a superfluid phase for weak disorder. The crit- 
ical exponents we estimate for this special multicritical 
point are different from those at the generic BG-SG tran- 
sition at incommensurate boson densities and agree with 
those for the pure three-dimensional XY model. This 
suggests that the latter and the tip of the lobe at weak 
disorder are within the same universality class. Rcnor- 
malization group and scaling arguments put forward in 
Q give strong support to this scenario. 
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